library(splm)
library(plm)
library(foreign)

setwd("../Replication Archive/") # set working directory

# load data stacked for lagsarlm() estimations 
data_stacked <- read.dta("Dataset_final_stacked.dta")
attach(data_stacked)

# load data for splm(), sphet() estimations 
data_splm <- read.dta("Dataset_final.dta")
attach(data_splm)
bfsDummies <- as.factor(bfs_nr) 
timeDummies <- as.factor(year)
data_splm.plm <- plm.data(data_splm,index=c("bfs_nr","year"))


bfsDummies <- as.factor(bfs_nr) 
timeDummies <- as.factor(year)

# Loading Connectivity Matrices
#------------------------------------------------#

# Time-invariant Structural Equivalence 1980
# Large Matrix
SEQ1980 <- read.dta("SEQ1980_Mat.dta")
w1980 <- as.matrix(SEQ1980)
weights_SEQ1980 <- mat2listw(w1980, row.names=NULL)
weights_SEQ1980 <- nb2listw(weights_SEQ1980$neighbours, style="W", zero.policy=TRUE)

# Short Matrix
s_weights1980 <- read.dta("SEQ1980_shortMat.dta")
s1980 <- as.matrix(s_weights1980)


# Planning Regions
# Large Matrix
weights_PR <- read.dta("PlanRegions_Mat.dta")
www_PR <- as.matrix(weights_PR)
weights_PR<- mat2listw(www_PR, row.names=NULL)
weights_PR <- nb2listw(weights_PR$neighbours, style="W", zero.policy=TRUE)

# Short Matrix
s_PR <- read.dta("PlanRegions_shortMat.dta")
sPR <- as.matrix(s_PR)


# Neighbours
# Large Matrix
weights_N <- read.dta("Neighbours_Mat.dta")
www_N <- as.matrix(weights_N)
weights_N <- mat2listw(www_N, row.names=NULL)
weights_N <- nb2listw(weights_N$neighbours, style="W", zero.policy=TRUE)

# Short Matrix
s_neigh <- read.dta("Neighbours_shortMat.dta")
sNeigh <- as.matrix(s_neigh)


# Hausman Test
#------------------------------------------------#

# Formulas for spatial Hausman Tests and later robustness checks using spml()
# (without unit dummies)

spml_nLDV <-   log_TotQuota ~ sh_Own   +
                sh_GreenParties         +
                sh_ConstrBus            +
                sqrt_spExpend_pKc       +
                sh_seFinance            +
                tax_r                   +
                log_popDensity_ob       +
                av_dPop5                +
                sqrt_ConstrSur_uz       

spml_LDV <-     log_TotQuota ~ Llog_TotQuota +
                sh_Own                  +
                sh_GreenParties         +
                sh_ConstrBus            +
                sqrt_spExpend_pKc       +
                sh_seFinance            +
                tax_r                   +
                log_popDensity_ob       +
                av_dPop5                +
                sqrt_ConstrSur_uz       

mod1.1 <- spgm(spml_nLDV, data=data_splm.plm, listw=mat2listw(s1980), lag = TRUE,
                moments = "fullweights", model = "random", spatial.error = TRUE)
mod1.2 <- spgm(spml_nLDV, data=data_splm.plm, listw=mat2listw(s1980), lag=TRUE,
                 model = "within", spatial.error = TRUE)
htest1 <- sphtest(x = mod1.1, x2 = mod1.2)
htest1

mod2.1 <- spgm(spml_LDV, data=data_splm.plm, listw=mat2listw(s1980), lag = TRUE,
             moments = "fullweights", model = "random", spatial.error = TRUE)
mod2.2 <- spgm(spml_LDV, data=data_splm.plm, listw=mat2listw(s1980), lag=TRUE,
             model = "within", spatial.error = TRUE)
htest2 <- sphtest(x = mod2.1, x2 = mod2.2)
htest2

# FE appropriate specification



# Table 5: Results
#------------------------------------------------#

# Formulas for main models estimated with lagsarlm()

nLDV       <- log_TotQuota ~ sh_Own   +
              sh_GreenParties         +
              sh_ConstrBus            +
              sqrt_spExpend_pKc       +
              sh_seFinance            +
              tax_r                   +
              log_popDensity_ob       +
              av_dPop5                +
              sqrt_ConstrSur_uz       +
              bfsDummies    

LDV        <- log_TotQuota ~ Llog_TotQuota +
              sh_Own                  +
              sh_GreenParties         +
              sh_ConstrBus            +
              sqrt_spExpend_pKc       +
              sh_seFinance            +
              tax_r                   +
              log_popDensity_ob       +
              av_dPop5                +
              sqrt_ConstrSur_uz       +
              bfsDummies   

LDV_noTaxR   <- log_TotQuota ~ Llog_TotQuota +
              sh_Own                  +
              sh_GreenParties         +
              sh_ConstrBus            +
              sqrt_spExpend_pKc       +
              sh_seFinance            +
              log_popDensity_ob       +
              av_dPop5                +
              sqrt_ConstrSur_uz       +
              bfsDummies 


########################################
# Column 1 & 2: Structural Equivalence 
#
# Structural Equivalence 1980
SEQ_nLDV <- lagsarlm(nLDV, data=data_stacked, weights_SEQ1980
                               ,method="eigen", tol.solve=1.0e-11, 
                               control=list(fdHess=TRUE))
summary(SEQ_nLDV)
bptest.sarlm(SEQ_nLDV)

SEQ_LDV   <- lagsarlm(LDV, data=data_stacked, weights_SEQ1980
                             ,method="eigen", tol.solve=1.0e-11, 
                             control=list(fdHess=TRUE))
summary(SEQ_LDV)
bptest.sarlm(SEQ_LDV)

# Robustness Check w/o tax rate, not shown in Table

SEQ_LDV_noTaxR   <- lagsarlm(LDV_noTaxR, data=data_stacked, weights_SEQ1980
                      ,method="eigen", tol.solve=1.0e-11, 
                      control=list(fdHess=TRUE))
summary(SEQ_LDV_noTaxR)
bptest.sarlm(SEQ_LDV_noTaxR)
# Results are very similar w/o tax rate

########################################
# Column 3: Planning Regions
#
# 
WPR_LDV <- lagsarlm(LDV, data=data_stacked, weights_PR
                                ,method="eigen", tol.solve=1.0e-10, 
                                control=list(fdHess=TRUE))
summary(WPR_LDV)
bptest.sarlm(WPR_LDV)


########################################
# Column 4: Neighbours
#
# 
WN_LDV <- lagsarlm(LDV, data=data_stacked, weights_N
                           ,method="eigen", tol.solve=1.0e-11, 
                           control=list(fdHess=TRUE))
summary(WN_LDV)
bptest.sarlm(WN_LDV)


# Online Supplementary Material
#------------------------------------------------#


# Table1: Robustness Checks I: Spatial Panel MLE & IV Estimation with GS2SLS 
#------------------------------------------------#

########################################
# Column 1 & 2: IV Estimation with GS2SLS 
#
sl<-function(x) slag(x, s1980)
ck<-function(x) crossprod(kronecker(x, diag(mat2listw(s1980))), x)

IVnoLDV <- plm(log_TotQuota ~   sl(log_TotQuota) +
# IVwLDV <- plm(log_TotQuota ~   sl(log_TotQuota) +
#                Llog_TotQuota           +                
                 sh_Own                  +
                 sh_GreenParties         +
                 sh_ConstrBus            +
                 sqrt_spExpend_pKc       +
                 sh_seFinance            +
                 tax_r                   +
                 log_popDensity_ob       +
                 av_dPop5                +
                 sqrt_ConstrSur_uz       |
#                Llog_TotQuota           +  
                 sh_Own                  +
                 sh_GreenParties         +
                 sh_ConstrBus            +
                 sqrt_spExpend_pKc       +
                 sh_seFinance            +
                 tax_r                   +
                 log_popDensity_ob       +
                 av_dPop5                +
                 sqrt_ConstrSur_uz       +
#                sl(Llog_TotQuota)       +
                 sl(sh_Own)              +
                 sl(sh_GreenParties)     +
                 sl(sh_ConstrBus)        +
                 sl(sh_seFinance)        +
                 sl(tax_r)               +
                 sl(log_popDensity_ob)   +
                 sl(av_dPop5)            +
                 sl(sqrt_ConstrSur_uz)             
               , data_splm.plm, model="within", effect="individual", inst.method="baltagi"); summary(IVnoLDV); fixef(IVnoLDV)
summary(IVnoLDV)
summary(IVwLDV)

########################################
# Column 3 & 4: Structural Equivalence using Spatial Panel MLE with spml()
#
Res_spml_nLDV <- spml(spml_nLDV, data=data_splm.plm, listw=mat2listw(s1980), model="w", effect="individual", lag=TRUE)
summary(Res_spml_nLDV)

Res_spml_LDV <- spml(spml_LDV, data=data_splm.plm, listw=mat2listw(s1980), model="w", effect="individual", lag=TRUE)
summary(Res_spml_LDV)

# Column 5: Spatial GLM with Random Effects

Res_spml_LDV_RE <- spgm(spml_LDV, data=data_splm.plm, listw=mat2listw(s1980), lag = TRUE, moments = "fullweights", model = "random", spatial.error = TRUE)
summary(Res_spml_LDV_RE) 


# Table 2: Robustness Checks II: Results with Period Fixed Effects
#------------------------------------------------#

########################################
# Column 1 & 2: Structural Equivalence w/ Unit & Period FEs using Spatial Panel MLE with spml()
#
spml_nLDV_TPFE <- spml(spml_nLDV, data=data_splm.plm, listw=mat2listw(s1980), model="w", effect="twoways", lag=TRUE)
summary(spml_nLDV_TPFE)

spml_LDV_TPFE <- spml(spml_LDV, data=data_splm.plm, listw=mat2listw(s1980), model="w", effect="twoways", lag=TRUE)
summary(spml_LDV_TPFE)

########################################
# Column 3: Regional Planning Groups w/ Unit & Period FEs  using Spatial Panel MLE with spml()
#
PR_spml_nLDV_TPFE <- spml(spml_nLDV, data=data_splm.plm, listw=mat2listw(sPR), model="w", effect="twoways", lag=TRUE)
summary(PR_spml_nLDV_TPFE)


########################################
# Column 4: Neighbourhood w/ Unit & Period FEs using Spatial Panel MLE with spml()
#
Neigh_spml_nLDV_TPFE <- spml(spml_nLDV, data=data_splm.plm, listw=mat2listw(sNeigh), model="w", effect="twoways", lag=TRUE)
summary(Neigh_spml_nLDV_TPFE)


